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Abstract 

Correlation analysis is convenient and frequently used tool for investigation of time 
series from complex systems. Recently new methods such as the multifractal detrended 
fluctuation analysis (MFDFA) and the wavelet transform modulus maximum method 
(WTMM) have been developed. By means of these methods (i) we can investigate 
long-range correlations in time series and (ii) we can calculate fractal spectra of these 
time series. But opposite to the classical tool for correlation analysis - the autocor- 
relation function, the newly developed tools are not applicable to all kinds of time 
series. The unappropriate application of MFDFA or WTMM leads to wrong results 
and conclusions. In this article we discuss the opportunities and risks connected to the 
application of the MFDFA method to time series from a random number generator and 
to experimentally measured time series (i) for accelerations of an agricultural tractor 
and (ii) for the heartbeat activity of Drosophila melanogaster. Our main goal is to 
emphasize on what can be done and what can not be done by the MFDFA as tool for 
investigation of time series. 
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Abstract 



KopejiaLi,HOHHHaT anajins e yflo6eH h hbcto HsnojiSBan HHCTpyMCHT sa nscjiefl- 
Baiie Ha BpeMeBH pe^OBe ot cjiojkhh CHCTeMH. B nocjiefliiHTe roflHHH 6Hxa pasBHTH 
HOBH MeTOflH 3a aHajiH3 Ha Kopejiai],HH KaTO Hanpniviep MyjiTHcJppaKTajiHHaT cJjjiyKTya- 
ii,HOHeH aHajiH3 (MFDFA) hjih meTOfl^vr na MaKCHMyMHTe na MOflyjiHTe na wavelet 
-TpaHc4)opMaLi,HHTa (WTMM). Hpe3 Te3H MCTOflH (i) nne MoaceM fla HSCJieflBaivie ko- 
pejian,HHTe c /i;'BjTL>r o6cer Bt>B bpcmcbh pe/i;oBe h (ii) fla HSHHCJiHBaMe (J)paKTajTiiHTe 
cneKTpH Ha BejiHHHHHTe, xapaKTepH3Hpan],H tcsh BpeMSBH peflOBe. Ho npoTHBHO na 
cjiynaa c KjiacHnecKHH HHCTpyMSHT 3a KopejiaLi,HOHeH aHajiH3 - aBTOKopejiai],HOHHaTa 
(J)yiiKn,HH, ropecnoMBHaTHTe MeTO/i,H iie ca npnjrojKHMH 3a BCsiKaKBH BpeMeBH pe/i,OBe, 
a 5e3pa36opHOTO hm npHJiarane MOJKe /i,a /i,OBe/i,e /lo rpeniHH pe3yjTTaTH h 3aKJTio^ieHH5i. 
B Ta3H CTaTHH HHe oScfcJKflaMe B'fc3MOJKHOCTHTe H pHCKOBeTe, CB'fcp3aHH c npHjiaraneTO 
Ha MFDFA MeTOfla kt^m BpeMeBH peflOBe ot e^HH renepaTop na cjiyHafiHH HHCJia h kt^m 
eKcnepHMeHTajiHO H3MepeHH BpeMeBH pe^OBe ot (i) ycKopennH na TpaxTop h (ii) ot 
CT.p/i;eHHaTa /i;eHHOCT na Drosophila melanogaster. HamaTa rjiaBna Ti;eji e fla HOflnep- 
TaeM KaKBO MOJKe h KaxBO ne Moace fla ce nanpaBH c MFDFA MeTOfla KaTO HHCTpyMenT 
3a H3CJieflBaHe na BpeMeBH peflOBe. 
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1 Introduction 



The investigation of correlation properties is one of the first tasks that a researcher of a time 
series has to perform. The classical tool for this is the autocorrelation function which can 
be calculated for stationary as well as for the nonstationary time series. The autocorrelation 
function reflects the behaviour of the time series. If for an example a time series is periodic 
its autocorrelation function is periodic too. For the most experimentally measured time 
series however the correlation function decays and in simple systems this decay in the most 
of the cases is exponential one. Opposite to this in numerous complex systems, which consist 
of many interacting subsystems, the correlations decay with power law and are called long- 
range because at large time scales the power law function is always large than the exponential 
function. 

The properties of the long-range correlations can be investigated by means of two recently 
developed methods: WTMM and MFDFA. WTMM is based on the properties of the wavelets 
[H El El E! and it leads to excellent results when the recorded time series are long enough. 
For shorter time series (usually below 10^ values) more convenient is the much simpler for 
implementation MFDFA method based on the scaling properties of the so called fluctuation 
function. Below we shall discuss what can be done and what can not be done by the MFDFA 
method when it is applied to different kinds of time series. We use two kinds of experimentally 
recorded time series (i) three short time series for the acceleration of an agricultural tractor 
working in different regimes and (ii) two time series for the heartbeat activity of Drosophila 
melanogaster. In addition we shall investigate a time series obtained by a random numbers 
generator in order to illustrate the monofractal properties of a time series. 

The paper is organized as follows. In the next section we discuss some multifractal 
quantities and how they can be calculated by the MFDFA method. In section 3 we give 
more information about the time series and analyze their autocorrelation functions. In 
section 4 we analyze the time series by means of the MFDFA method. Several concluding 
remarks are summarized in the last section. 

2 Generalized dimensions, fractal spectra, Holder and 
Hurst exponents 

Let us have a set S which lies in an A^— dimensional Cartesian space covered by a grid of 
A^-dimensional boxes of edge length e. Let e be small and we need A^*(e) boxes to cover our 
set. The box-counting dimension of the set is 



Dq is not an integer for some sets. These sets are called fractals. Dq is a member of the 
spectrum of the Dg dimensions 



(1) 




(2) 




where g is a continuous index and I{q,e) 



N*{e) 



(3) 
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is a sum over the N* boxes of size e which cover the set. /i^ is the natural measure i.e. it is a 
measure of the frequency with which typical orbit visit various boxes covering the attractor 
for the limit case when the length of the orbit goes to infinity (in addition the frequencies 
have to be the same for all initial conditions in the basin of attraction of the attractor except 
for a set with Lebesque measure 0). Thus for fi^ we have 

4 Uh = lim — , 

where ^ is the time the orbit originating from xq spends in the cube in the time interval 
< t < T. From ^ by means of the L'Hospital rule we easy can obtain 

(5) A = lim.^^-^^)^^^"^^ 



-0 Ine 

Di is called also information dimension. In general Dq > Di > D2 > . . .. If -Dg varies with 
q the measure, associated with Dq is called multifractal measure. The measure can have 
different behavior in each cube covering the set regardless of how small is its size e. Thus 
the dimension is not enough for a characterization of the measure and we need additional 
characteristic quantities. One such quantity is the coarse-grained Holder exponent for each 
box 

, , \nu(box) 

/X can have irregular behavior for different boxes. When the box size decrease and /i has 
statistically the same irregular behavior the measure is called self-similar jH] . For large class 
of such measures a lies between Omm > and amax < 00. For given e we can count the 
number of the boxes with value of a in small interval around a prescribed value. Thus 
we can define the coarse-grained spectrum 

(7) Ma) oc 

From (|7j) we obtain A^^^ (a) ^ e-^("), i.e. the relationship which has a form analogous to this 
one for the dimension above. Thus /(a) can be treated as fractal dimension of the subsets 
Sa of S having coarse grained Holder exponent equal to a. If /(a) is not constant then for 
small e the set S is constructed by subsets Sa of different dimensions /(a). Then the set S 
is called multifractal. 

The multifractals can be characterized by different spectra. Below we shall obtain rela- 
tionships for these spectra for the case e — *• 0. We remember that the set 5* be covered by 
grid of boxes of unit size e and n is the probabihty measure on 5* {fi{S) = 1). Let /i(cfc) = fik 
where Ck denotes again the k—th cube. We have assigned a singularity measure ak to each 
cube 

(8) f^k = e"'= 

For small e we can make continuous approximation for the number of cubes for which ak is 
between a and a + da, i.e., we can denote this number as 

(9) p{a)e-^^"^da 
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By substitution of (jH)) in the relationship for /(g, e) and after a transition from a sum over 
the boxes to an integration over the a we obtain 



fc=i 

(10) = y exp {[/(a*) - ga*] ln(l/e)} 

For small e ln(l/e) is large and the main contribution to the above integral is from the 
neighborhood of the maximum value of the /(a*) — qo* . Let f{a*) be smooth and the 
maximum is located at a* = a{q) given by 

(11) ~ \a'=a{q)= ^ ^ l<^*=«= Q 

(12) :7r^[/(«*)-^«1 ■ "^'^ ' 



Now we take the Taylor series representation of the function F{a*, q) = /(a*) — qa* around 
the point a* = a{q) and substitute it in (fTUj). The result is 

I{q,e)=exp{[f{a{q))-qa] ln(l/e)} J rfa*p(a*)e-(^/2)/" K,))(a*-a(,))2 

(13) ^exp{[/(a(g))-ga]ln(l/e)} 
Introducing in ^ we obtain 

(14) ^'^^ ^^""^^^ " ^^"^^^^^ 
From dni) 

(15) = „(,) = I 
Then 

(16) r(g) = (g - 1)D, ^ D, - ^ 



q-1 



From (HD) 



(17) /(«(g)) = g— - (^ - 1)^. = - ^(^) 



dr dr 

Thus for each g from (|TT)|l and (fTTjl give a(g) and /(a) thus parametrically specifying the 
function /(a). The Hurst exponent sometimes is associated with the coarse-grained Holder 
exponent for e — > 0. This can be done in the following way. For the cases in which the 
following relationship holds 

(18) r(g) = qh{q) - 1 
we obtain 

dr , , , dh 

(19) = 

and 

(20) /(a) = ga - r(g) = q[a - h{q)] + 1 
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But in which cases holds fITHj)? Let us consider stationary time series {xk}, k = 1, . . . , N and 
let us use appropriate transformations in order to make all values positive and to normalize 
the time series, i.e., 

N 

(21) Xk>0, E^fc = l 

k=l 

In this case we can associate the time series with some probabilities and we shall use this 
fact to derive ()18j) . Let us construct the profile function for our normalized time series 

n 

(22) Y^ = J2{xk-{x)) 

k=l 

where (x) is the mean of the time series. Now we divide the time series into Ns segments of 
length s and let for simplicity Ns = N/ s is an integer. The sum 

us 

(23) Yius)-Y{iu-l)s)= ix,-{x)) 

k={u-l)s+l 

for the segment u is identical to the box probability Ps{i^) which is the main building block 
of the corresponding partition sum Zq{s) and scales with the Tq spectrum 

N/s 

(24) Zqis) = Y.\ps{u) |''ocs-('^) 

u=l 

From other side for our time series {xk} we can build a scaling function which contains the 
local Hurst exponent. This is the sum 

(25) ^Ein^^)-n(^-l)Hn ^ 



From here we obtain 



2N/s 



(26) -^E |r(z..)-r((z.-i).r(xs«^(«) 

and finally 

N/s 

(27) E I F(z/s) - r((z/ - l)s l^oc s'?'^^'?)-^ 

A comparison of ^ and (j2I|) leads to (IT^ . 

The implementation of the MFDFA method j5j follows the steps below. The starting 
point are our time series {xk} of finite length N. First of all we have to determine the profile 
function. Here we have to mention the following. There are two possibilities. First of all we 
can use profile function as in (j22I), i-e., 

(28) Y{z) = ±{x,-{x)) 

k=l 

The final result of our analysis will be to obtain some values of the local Hurst exponent h. 
If we use ()28|1 as profile function we can determine only positive Hurst exponents which are 
not quite close to 0. In all other cases we have to use as profile 

(29) Y{{) = ±{Y{k)-{Y)) 

k=l 



where Y{k) comes from ()28|) . Thus we shall obtain value of local Hurst exponent which is 
larger than the value of true exponent i.e. 



(30) hiq) = 1 + h{q) 

This is the final difference between the two branches of the method. We shall present here 
the first branch, which has the same steps as the second with the only difference that we have 
everywhere to change Y with Y if we want to obtain the algorithm for the second branch. 
After the calculations of the profile function we have to divide it into segments of length s. 
As N/s in general is not an integer some part at the end of the time series will remain out 
of the Ns = int{N/s) segments. In order to incorporate the influence of this part into the 
analysis we divide the time series from the end to the beginning again to Ng segments of 
length s. Thus we obtain 2Ns segments and for each segment calculate the local polynomial 
trend yu{p), P = 1, . . . ,s. The trend can be a polynomial of order m. Thus we perform 
MFDFA of order m and the coefficients of the polynom are obtained by a least-square fit of 
the corresponding segment. After the trend calculation we determine the variance 

(31) F\u, s) = - ±{Y[{u - l)s + z] - yA^)r 

^ 1=1 

for the segments u = 1,2, . . . , Ng. For the segments u = Ns + l,Ns + 2, . . . , 2Ns the variance 
is 

(32) F\u, s) = - i:{Y[N - (z/ - N,)s + z] - y^iz)}' 

1=1 

The next step is an averaging over the segments for obtaining the g-th order fluctuation 
function [q is in general a real number). For g = the fluctuation function is 

(33) Fo{s)=expl^j:Hn^,s)]\ 



and for g 7^ 

(34) F,is) 



For large class of time series the fluctuation function has a power law behavior and from this 
power law for large s we can determine the local Hurst exponent 

(35) F,(s) cx s^(^) 



3 The time series and their autocorrelation functions 

In order to achieve our goals we have selected time series from three different systems: an 
agricultural tractor, time series from the heartbeat activity of the classical object of Genetics 
-Drosophila melanosgaster and a time series from a random number generator in order to 
illustrate the case of absence of correlations as well as the monofractal properties of some 
time series. The investigated time series are shown in Fig. 1. We investigate two kinds 
of experimental time series for acceleration of agricultural tractors. The time series from 
panels (a) and (b) of Fig. 1 are from experimental modeling of the bouncing phenomenon. 
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It arises when a machine travels on a rough road and consists of large-amphtude oscillations 
which can lead to injuries of the driver, to a lift of the machine (some times more than 
half-a-meter over the surface) and to damage of some of its elements. The time series for 
this phenomenon are recorded when the tractor moved on artificial rough road, consisting 
of 12 small kinks prepared on am asphalted surface. For this case we have several large 
amplitude oscillations when the tractor moves on the rough part of the road and in some 
time series we have also the noise background when the tractor moves on the smooth part 
of the road. Thus we can discuss these time series as time series of noise with superimposed 
periodic trend. The time series from panel (c) of Fig. 1 are recorded when the tractor has 
a construction exhibiting impact properties. This case has been chosen because the impact 
systems are frequent source of chaotic vibrations [ZI. The tractor has been converted 
to an impact system by adding vibrating subsoiler to it. The soil cutting chiesel breaks a 
hard soil layer located at the depth between 30 and 50 cm from the field surface and it is 
oscillated in order to reduce the draft force and to improve the water infiltration into the 
soil. The oscillation is realized by a hydravlic power cylinder with two lift arms. This work 
mode is called ground-penetrating mode. 

The time series in panels (d) and (e) of Fig. 1 are records of the heartbeat activity of 
Drosophila melanogaster These flies are provided by Bloomington Drosophila Stock Center, 
U.S.A. The time series are ECGs (electrocardiogramms) for the first generation (the kids) 
obtained by male Dopa decarboxilase (Ddc) mutant (FBgn 0000422) located in chromosome 
2, locus 37C1, crossed with female shibire (shi) (FBgn 0003392) located in chromosome 1, 
locus 13F7-12. Ddc codes for an enzyme necessary for the synthesis of four neurotransmitters: 
norepinephrine, dopamine, octopamine, serotonin, related to learning and memory. The 
shibire (shi) mutants cause paralysis at high temperature and eliminates the effect of the 
neurotransmitters on the heart. The ECGs have been recorded optically at a stage PI 
(white puparium) of a Drosophila development when it is both immobile and transparent 
and the dorsal vessel is easily viewed. The object was placed on a glass slide in a drop 
of distilled water under a microscope (magnification 350 x). Fluctuation in light intensity 
due to movement of the dorsal vessel tissue was captured by photocells fitted to the one 
eyepiece of the microscope. The captured analogue signal was then digitized at 1 kHz 
sampling rate by data acquisition card and Lab VIEW data capturing software supplied 
by National Instruments. Finally the time series from the random number generator are 
obtained by means of the generator ran2 from JUj. In Fig. 2 are presented the autocorrelation 
functions for the investigated time series. The large amplitude oscillations in panels (a) and 
(b) reflect the bouncing phenomenon. The impact oscillations are clearly reflected by the 
almost periodic behaviour of the autocorrelation function in panel (c). The slowly decaying 
autocorrelations for the ECGs of Drosophila show some degree of anticorrelation at large 
n. Finally as it can be expected the time series from the random number generator are 
uncorrelated. 

4 Fractal analysis of the time series 

Excluding the time series from the random number generator which are uncorrelated, the 
other time series exhibit some degree of long range correlations. Below we shall discuss the 
question how far these correlations can be investigated by means of the modern methods 
such as the MFDFA. In order to answer this question we have to look at the fluctuation 
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functions connected to the time series. We shall calculate the fluctuation functions directly 
on the basis of the time series for the cases of tractor accelerations and for the time series 
from the random number generator. For the time series from Drosophila we first shall 
construct the time series for the intermaxima intervals (the time series for the time between 
two consequent maxima of the time series) and for the intermaxima time series we shall 
calculate the fluctuation functions and eventually the characteristic fractal quantities. The 
results are shown in Fig. 3. For the time series for the acceleration of the tractor in bouncing 
regime we can not apply the MFDFA method. The reason is that these time series (shown 
in panels (a) and (b) in Fig. 1) are too short and we can not divide them in large number of 
segments in order to have enough statistics for correct calculation of the fluctuation function. 
In order to apply the MFDFA the time series must have at least two thousand values. The 
minimum number of values - 4000 for which the fluctuation function is calculated here are 
the time series for the acceleration of the tractor in impact regime. The fluctuation function 
is not a straight line on a log-log scale - it has a kink as a reflection of periodicity of the time 
series. Thus these time series do not show scaling behaviour and we can not calculate the 
fractal spectra for them. On the other side we confirm the observation [Hj that the periodic 
behaviour influences the fluctuation function. On panels (b) and (c) of Fig. 3 are presented 
the fluctuation functions for the intermaxima intervals of Drosophila. On panel (b) the 
fluctuation function is close to a straight line, i.e., the scaling properties can be assumed 
and the multifractal spectra can be calculated. Not so clear is the situation in panel (c). 
We shall carry further the calculations for this case in order to see how the not very large 
deviations from the scaling behaviour influence the form of the fractal spectra. The best 
scaling properties are exhibited by the time series from the random number generator - 
panel (d). This can be expected as these time series are specially selected for illustration of 
monofractal behaviour. Thus for calculation of the fractal characteristics from the initial 6 
time series the requirements for length and presence of scaling properties have eliminated 3 
time series. The h{q) spectra for the remaining three time series are presented in Fig. 4. 
The expected behaviour of such kind of spectrum is presented in panel (a) where h decreases 
with increasing q. This multifractal behaviour reflects the scaling properties of the time 
series as shown in panel (b) of Fig. 3. Multifractal behaviour is exhibited also by the second 
time series for the intermaxima intervals of Drosophila melanogaster shown in panel (b) of 
Fig. 4. As it can be seen from panel (c) of Fig. 4 for the time series from the random 
number generator h is a constant for all q which is characteristic feature for the monofractal 
behaviour. 

Finally we calculate the /(«) spectra. For the time series from the random number 
generator this spectrum is a point and it is not presented in Fig. 5. Panel (a) shows the 
spectrum for the time series with scaling properties shown in panel (b) of Fig. 3. The 
expected form for the /(a) spectrum is parabolic one. This form is observed on panel 
(a) of Fig. 5. The maximum of the /(a) spectrum shows at which a is positioned the 
most statistically significant part of the time series, i.e., the subsets with maximum fractal 
dimension among all subsets of the time series. In our case fmax = 1 at a = 0.652. These 
and other characteristics of the f{a) spectrum (such as its width at some value k fmax j 
< k < 1) can be used as tools for classification of time series with multifractal properties. 
The consequences of not fully scaling properties of the time series from panel (c) of Fig. 3 are 
observed in panel (b) of Fig. 5 where the parabolic form of the spectrum is slightly deformed. 
As such a form is different from the form for which we can make quantitative conclusions 
we can use such spectra only for qualitative conclusions and classifications. From the point 
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of view of genetics the different spectra correspond to the different genetic characteristics 
which the flies from the flrst generation obtained from the parents with different genetic 
heart defects. Thus only on the basis of the ECGs we can conclude that the two investigated 
Drosophila flies are affected in different way by the genetic defects of their parents. Thus the 
heartbeat dynamics of Drosophila is connected to its genetics. This conclusion opens the way 
for further intensive research on the relation between (i) the biomechanical properties of the 
heart of simple and more complex animals and humans and (ii) their genetic speciflcations. 

5 Concluding remarks 

In this paper we demonstrate the possibilities and risks when the multifractal detrended 
fluctuation analysis is applied to experimentally obtained time series. We have started with 
six time series from which only one was carefully selected in order to demonstrate monofractal 
and irregular behaviour. Then step by step we have demonstrated how the requirements of 
the method stop the investigation of different time series at different points as follows. The 
requirement for length was not satisfled by the time series for the acceleration of the tractor 
in bouncing regime. The investigation of the time series for the tractor acceleration in impact 
regime was stopped at the point where the requirement for scaling of the fluctuation function 
was imposed. At this point we allowed us to make a relative crude approximation for scaling 
properties of one of the time series for the heartbeat activity of Drosophila. The above crude 
approximation has the consequences in the multifractal spectra and especially in the /(«) 
spectra. The monofractality of time series was clearly illustrated by the h{q) spectrum of 
the time series from the random number generator. Finally the f{oi) spectra have shown 
that the prescribed form of the spectrum for the multifractal behaviour can be assigned to 
only one of the six time series. 

In conclusion we warn the researcher to be very careful when he or she applies the MFDFA 
to time series from real systems. Only small number of them are appropriate for such 
analysis. But when this analysis can be performed it can supply us with much quantitative 
and qualitative information about the relations between (i) the observed dynamics of the 
investigated complex system and (ii) the usually not visible processes, responsible for this 
dynamics. 
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Figure captions 

Figure 1: 

The time series. Panels (a), (b), (c): time series for the acceleration of an agricultural 
tractor. Horizontal axis : time (unit for time is 8 ms for the panels (a) and (b) and 
1 ms for the panel (c)). Vertical axis : accelerations (unit m/s^). Panels (d) and (e): 
characteristic time series for the heartbeat activity (ECG) of Drosophila melanogaster. 
Unit for time here is 1 ms. Panel (f): time series from the random number generator 
ran2 ^U]. This generator generates pseudorandom numbers in the interval [0,1]. 

Figure 2: 

Autocorrelation functions for the time series shown in Fig.l. Panels (a), (b): auto- 
correlations for the bouncing tractor. Panel (c): autocorrelations for a tractor work- 
ing in impact regime. Panels (d), (e): autocorrelations for the heartbeat activity of 
Drosophila melanogaster. Panel (f): time series of pseudorandom numbers obtained 
by ran2. 

Figures : 

Fluctuation functions ()S4|1 . For all panels from bottom to the top the curves (marked 
with circles) are obtained for q = 2,4,6,8. Solid line denote the r.m.s. power law 
fits of the corresponding curves. Panel (a): fluctuation functions for the time series of 
agricultural tractor in impact regime. Panels (b) and (c): fluctuation functions for in- 
termaxima time series obtained by the time series of ECGs of Drosophila melanogaster. 
Panel (d): fluctuation function for the time series obtained by the generator ran2. 

Figure 4: 

Local Hurst exponent spectra h{q). Panel (a): h{q) spectrum for the intermaxima 
time series obtained by the ECG time series of Drosophila shown in panel (d) of Fig.l. 
Panle (b): h{q) spectrum for the intermaxima time series obtained by the ECG time 
series of Drosophila shown in panel (e) of Fig. 1. Panel (c); h{q) spectrum for the time 
series obtained by the random number generator ran2. 

Figure 5: 

/(a) spectra for the intermaxima time series obtained by the ECG of Drosophila 
melanogaster. Panel (a): /(a) spectrum for the intermaxima time series oftained by 
the ECG of panel (d) of Fig. 1. Panel (b): /(a) spectrum for the intermaxima time 
series obtained by the ECG of panel (e) of Fig. 1. 
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Fig. 4 
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